Electronic polarons and bipolarons in Fe-based superconductors: a pairing mechanism 
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Superconductivity is a fascinating example of 
how "more is different" [1]. It is due to electrons 
binding into bosonic Cooper pairs, which exhibit 
coherent behavior across a macroscopic sample. 
Finding the mechanism responsible for this bind- 
ing is one of the more difficult tasks of condensed 
matter physics. For conventional superconduc- 
tors the solution was given by the BCS the- 
ory [2] as being due to exchange of phonons. For 
the cuprate high-T c superconductors a widely- 
accepted explanation is still missing despite in- 
tense effort. The recently discovered Fe-based 
high-T c superconductors pose now a new chal- 
lenge [3, 4, 5, 6, 7]. We present here a quantum 
mechanical theory for pnictides describing the in- 
fluence of the large electronic polarizability of the 
heavy anions. We demonstrate that its inclusion 
results in electronic polarons as the low-energy 
quasi-particles and also unveils a pairing mecha- 
nism for these electronic polarons. 

At first sight, one may expect the pairing mechanism 
of pnictides to be related to that of the cuprates, given 
the somewhat similar layered structures. However, there 
is an essential difference between the Cu02 layer where 
the important physics takes place in a cuprate, and its 
counterpart, the FeAs layer of a pnictide: while the for- 
mer is two-dimensional (2D), with all Cu and atoms in 
the same layer, the latter is not. Instead, the layer host- 
ing the Fe atoms (which are arranged on a simple square 
lattice) is sandwiched between two layers which share 
equally the As atoms, as sketched in Fig. 1(a). Each 
Fe has 4 nearest neighbor (nn) As atoms at a distance 
R = 2.4A, arranged in a somewhat distorted tetrahedron, 
two in the upper and two in the lower layer. 

The different geometry has important consequences. 
In the cuprates, the states near the Fermi energy consist 
of strongly hybridized Cu 3d and O 2p orbitals. In con- 
trast, hardly any hybridization appears between the Fe 
3d and As 4p orbitals, with the former giving essentially 
all the contribution to the low-energy states of the pnic- 
tides [8]. This lack of hybridization is due not only to 
the lattice structure but also to the substantial spread in 
space of ip orbitals, as compared to 2p orbitals. 

As a result, one might assume that the simplest Hamil- 
tonian describing the low-energy physics of pnictides is a 
Hubbard Hamiltonian for the Fe 3d electrons, namely: 
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FIG. 1: (a) 3D sketch of the idealized FeAs layer. The lattice 
constant a — 2.8A and angle 8 (see text) are indicated; (b) 
Top view of the same Fe-As layer. Several sites are indexed. 
The Fe and As atoms are at different "depths". 



where c] a creates an electron on the Fe site i with spin 

a and hi„ = c^ ia Ci a . For simplicity, here we only count 
"doping" electrons, i.e. extra charges on top of the 3d 6 
configuration of the Fe in the undoped compound. Of 
course, a multiple-band Hamiltonian can also be used for 
the proper description of the various 3d orbitals, how- 
ever the essential physics we want to discuss is already 
captured within this simpler starting point. The hopping 
integral is certainly finite for nn Fe sites, with ti^ +x = t. 
We will also consider the effects of including 2 nd nn hop- 
ping, with U,i +x+y — t' [the indexing of various sites is 
indicated in Fig. 1(b)]. The hopping to As sites is ne- 
glected, given the small hybridization. 

Of course, Hamiltonian (1) is the most common start- 
ing point for cuprates (where it does include the O contri- 
bution, as c\ a are operators for Zhang- Rice singlets [9]). 
In contrast, using (1) or a similar starting point either 
in the strong-coupling [10, 11, 12, 13, 14, 15] or weak- 
coupling limit [16, 17, 18, 19] as the low-energy Hamilto- 
nian for pnictides implies that the As anions play no role 
in their physics. 

However, the As atoms are not irrelevant. As pointed 
out recently in Ref. 8, these are big, highly-polarizable 
ions which are strongly influenced by the extra charges 



in their vicinity. One expects each such charge to be sur- 
rounded by polarized As ions, giving rise to electronic 
polarons. This results in a strong screening of the on- 
site Coulomb repulsion, suggesting that these materials 
are not in the large U limit of a Mott-Hubbard insu- 
lator. It also results in a strong nn attraction, which 
may be the key component in the pairing mechanism. 
The arguments of Ref. 8 are based on semi-classical esti- 
mates. Here we propose a quantum model which allows 
us to investigate this phenomenology away from the lin- 
ear regime, and to also study dynamic properties of these 
polarons and bound bipolarons. 

The polarization of the As 3 " ions by the electric fields 
of the external charges is due to matrix elements cou- 
pling their filled Ap to (primarily) their empty 5s orbitals 
(more details are in the supplementary material). This 
is equivalent to a hole being virtually excited from the 5s 
into the 4p orbitals. We therefore describe each As by 4 
operators: s\,p\. a ,v\ <jiV\ <r> which create a hole of spin 
a in the respective orbital. The ground-state of an As 
ion is sjsj |0) and the p-orbitals lie at an energy f2 above 
it, so that the unperturbed As ions are described by: 
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where i is the location of the As ion [see Fig. 1(b)] 
X — x,y,z and a is the spin of the hole. Hopping between 
As ions is very small, given the large nn distance between 
them of about 4A, and we ignore it. Indeed, LDA calcu- 
lations find narrow Ap bands [8] , justifying this atomistic 
description of the As ions. 

We make two approximations regarding the interac- 
tions: (i) only As atoms nn to an Fe hosting a charge are 
polarized by its electric field, and (ii) we ignore dipole- 
dipole interactions between As clouds. (Estimates of the 
corresponding corrections are given in the supplementary 
material). Our interaction Hamiltonian is, then: 
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where the sum is over all unit cells i in the lattice. This 
interaction describes the 5s <-> Ap transitions resulting 
in the polarization of the As along the appropriate Fe- 
As direction (see Fig. 1) if doping charges, counted by 
hi = rijf + hii, are on any Fe^ site. 

Our model Hamiltonian H — Hp c + Has + Hint is 
characterized by 5 energy scales. The hopping integral 
t = 0.25eV and the 5s — Ap energy difference f2 = 6eV 
are estimated from LDA results (for details on determin- 
ing all the parameters, sec the supplementary material). 




FIG. 2: Sketch of a single electronic polaron and its nn hop- 
ping. The electric field created by the extra electron residing 
on the Fe excites holes of nn As into the Ap orbital pointing 
towards the Fe (more distant As atoms stay in the 5s ground 
state). The hopping integral is renormalized by the overlap 
of the polarization clouds in the initial and final states. 
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t/2, the latter being the 
or- 
bitals which contribute most near the Fermi level. The 
interaction energy g = 2.5eV is extracted from the As 3- 
polarizability a p = 10A 3 , and the Hubbard onsite repul- 
sion Uh is left as a free parameter. We present results 
for a wide range of these parameters, proving that our 
results are only weakly sensitive to their precise values. 

We begin with a study of the low-energy spectrum of 
the electronic polaron [sketched in Fig. 2] created when 
there is one doping electron in the system. We use first 
order perturbation theory in the hopping, which is the 
smallest energy (a quantitative justification and full de- 
tails of the calculations are provided in the supplemen- 
tary material). The polaron eigenenergies are: 
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The first term is the interaction energy between the elec- 
tron and the induced As dipole moments. For g <C ^ 
it equals — A[a p E 2 /2] «-8.2eV, E being the electric field 
at the As sites. This linear expression was used in Ref. 
8. Since in fact g/il ~ 0.4, non-linear effects become 
important and reduce it slightly to ss — 7.1eV. e c g(k) = 
—2t c g \cos(k x a) + cos(k y a)] — At' oS cos(k x a) cos(k y a) is a 
dispersion identical to that of a free particle, but with 
renormalized hopping integrals due to the overlap of the 
As clouds as the polaron moves between different sites 
(see Fig. 2). Expressions for t e g and t' eS are given in 
Eqs. (23), (24) in the supplementary material. We plot 
their values in Fig. 3, for a wide range of f2 and a p values. 

For a fixed a pi a smaller £1 implies a larger g/fl ra- 
tio [Eq. (14) in supplementary material], i.e. a stronger 
effective coupling. This explains the smaller t c g, t' oS as 
f2 decreases to be typical polaronic physics: stronger in- 
teractions lead to heavier polarons [20]. However, note 
that these ratios are not exponential in the coupling, like 
for normal lattice polarons, but rather weakly dependent 
on the various parameters. This is due to non- linear ef- 
fects, namely the existence of a maximum allowed value 
of the induced dipole moment. Both t c s/t and t' cf{ /t' are 
around 50% for a wide range of values centered on our 
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FIG. 3: (a) t e s/t and (b) t' eS /t' vs. Q, for a polarizability 
a p = 7, 10 and 12 A 3 . The dots show the values used here. 



parameters, suggesting a light polaron with an effective 
mass roughly twice as big as the free band carrier. 

Such renormalization of the band-structure compared 
to that predicted by LDA should be detected by ARPES. 
Interestingly, a recent measurement [21] showed reason- 
able fits to LDA upon rescaling the bandwidth by a factor 
of 2.2, similar to our typical values (the material stud- 
ied has P instead of As, and these smaller anions have a 
smaller polarizability a p ss 7 A 3 ). To first order in per- 
turbation theory, we also predict a quasiparticlc weight 
which is independent of k, with Z — 0.38 for our param- 
eters. The remaining spectral weight should be observed 
at much high energies typical of the excited As polariza- 



tion clouds, of the order 
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higher (more details arc in the supplementary material). 

We now study a system with 2 doping electrons, also 
within first order perturbation in the hopping. We fo- 
cus on singlet solutions (eigenstates are either singlets 
or triplets of the two electrons). Also, we first analyze 
results for t' = 0, then investigate the role of a finite t'. 

If there were no interactions between the two polarons, 
then for a total momentum k one would get a continuum 
of eigenstates with energies Ep (k — q) + Ep(q) for all q in 
the Brillouin zone, with the ground-state at k = q = 0. 
However, there are short range interactions between the 
polarons. If the charges are 2 nd nn or closer to each other, 
some As ions are simultaneously nn to both charges and 
acquire a dipole moment different from that of As ions 
which interact with a single charge. This leads to differ- 
ent energies Uq, U\ and Ui for configurations with the 
two charges being, respectively, on-site, nn and 2 nd nn 
to each other. We sketch below the derivation of Uq and 
U\ . The effective hopping between these configurations is 
also different from t e s, because of the overlap with these 
differently polarized clouds. These interaction energies 
and some of the effective hoppings are summarized in 
Fig. 4 (expressions and plots for all these quantities are 
provided in the supplementary material). 

An on-site bipolaron also has 4 polarization clouds, 




FIG. 4: Sketches of on-site, 1 st , 2 nd and 3 rd nn bipolarons. 
The first three configurations have interaction energies Uo, Ui 
and U2, respectively. Several of the special effective hopping 
integrals are also indicated. 



however effectively g — > 2g in Eq. (3), as the on-site 
charge is doubled, so its electrostatic energy can be imme- 
diately obtained from Eq. (4) to be Ebp$ = Uh + A(H — 
^fl 2 + 16g 2 ). The energy of two static distant polarons 
is E B p,oc — 8(fl — \JVL 2 + 16g 2 ) , therefore Uq = E B p,o - 
v/ft 2 + IQg 2 +Q- 2^fl 2 + Ac 
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Uh always. This is obvious in the linear regime: because 
the electric field is twice as large, the energy of the bipo- 
laron is 4 times that of a polaron, i.e. twice as large as 
of two free polarons. This leads to a significant screening 
of Uh, of several eV, as shown in Fig. 5(a) (these values 
are smaller than in Ref. 8 because of non-linear effects) . 
Such a strong on-site attraction is well-known to arise for 
polarons in general [22]. 

A nn bipolaron has 6 polarization clouds (see Fig. 4), 
the two central ones being larger than the regular polaron 
clouds. In the linear regime, the energy of each central 
cloud is —a p (Ei + E2) 2 /2, so the energy difference with 
respect to two regular clouds is — a p E\ ■ E 2 < if the 
Fe-As-Fe angle is less than 90°, as is here. This means 
an attractive nn interaction, as indeed seen for typical 
parameters in Fig. 5(b). However, because of non-linear 
effects, this interaction becomes repulsive at strong cou- 
pling, as shown in Fig. 5(d). Interestingly, this figure 
reveals that these materials are close to optimal, with U\ 
near its minimum, for typical values of the parameters 
(see also Fig. 3 in the supplementary material). 

This nn attraction is not typical polaron physics. On a 
2D lattice like in cuprates, the interaction between elec- 
tronic polarons is repulsive because the Cu-O-Cu angle 
is 180° (the same holds true for polarons coupled to a 
breathing-mode phonon of the O atoms). For a Holstein 
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FIG. 5: (a) Renormalization of on-site interaction, Uo — Uh', 
(b) nn energy U\ and (c) 2 nd nn energy U2 vs. Q, for various 
polarizabilities. (d) Ui vs. g/Q, when Q. — 4, 6, 8eV. The dots 
show our typical values. 
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FIG. 7: Dispersion of the two bound bipolaron states along 
high-symmetry axes in the Brillouin zone, for (a) t' = and 
(b) t' — — 1/2. The two-polaron continuum is also shown. 
Parameters are U H = 10 eV, a p = 10 A 3 , fi = 6 eV (similar 
results are found for all a p = 7 - 12A 3 , fi = 4 - 8 eV). The 
symmetry of the ground state changes from s to d if t' 7^ 0. 



bipolaron, there is a weak nn attraction, but it is due to 
spin exchange and is 0(t 2 ), not O(t ) like here [22]. Its 
non-2D geometry is the essential ingredient in bringing 
about this strong nn attraction, for the Fe-based super- 
conductors. This mechanism is absent in cuprates. 

First order perturbation in hopping mixes these static 
configurations. Following the calculation described in 
the supplementary material, we obtain the eigenergies 
EBp(k) for a total momentum k. In Fig. 6(a) we plot 
k = eigcnenergies for all bound bipolaron states, i.e. 
which have energies below the two-polaron continuum 
starting at — 8i c ff- We find 3 such states. The GS first 
shows linear dependence on Uh an d then flattens out. 
Its wavefunction has s-type symmetry (unchanged sign 
upon 90° rotation). There is a second s-state at small 
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FIG. 6: (a) Eigenstates below the two-polaron continuum, 
AEbp = Ebp(0) + 8t a B vs. Uh- The symmetry of the three 
bound eigenstates is labeled in the inset, (b) Probability for 
on-site, 1 st , 2 nd and 3 rd nn separation in the ground state, 
vs. Uh- The dashed lines show the same quantities for the 
d-state. We use t' = 0,a p = lOA 3 , f2 = 6eV (similar results 
are found for all a p = 7 - 12A 3 , fi = 4 - 8 eV). 



Uh, which then joins the continuum. The 3 rd bound 
bipolaron state has an energy independent of Uh, and is 
d-type (wavefunction changes sign upon 90° rotation). 

The nature of these bound states is revealed in Fig. 
6(b). For low Uh values, the GS primarily consists of 
an on-site bipolaron, with hardly any contribution from 
nn or more distant configuration. This explains why its 
energy here scales with Uh (more precisely with Uo). 
When Uo « there is a fast crossover to a state dom- 
inated by the nn bipolaron configurations (a combined 
total of 90% probability). This is expected, since U\ < 
irrespective of Uh, favoring such a pair when Uo becomes 
repulsive. The onsite contribution is now exponentially 
small. This explains the weak Uh dependence here, as 
coming from virtual hopping to the on-site configuration. 
The dashed lines show the contributions for the d-type 
state. As expected, it is dominated by nn bipolarons. It 
has a zero on-site probability, consistent with its symme- 
try and explaining the lack of dependence on Uh- The 
2 nd nn contribution is also zero (U2 > as well). There 
are small 3 rd and more distant bipolaron contributions. 

The unscreened Hubbard repulsion is very large, Uh ~ 
10 — 20 eV. For these values there is hardly any depen- 
dence on Uh, so its precise value is of little importance. 
We use Uh = 10 eV from now, and ask how mobile are 
these bound, predominantly nn, bipolaron pairs. Their 
dispersion in the Brillouin zone in shown in Fig. 7(a), 
where we also show the two-polaron continuum (eigen- 
states above the continuum are not shown). The s-pair 
has a bandwidth Ecsi^, n) — Ecs(®, 0) ~ 0.3 eV, imply- 
ing a bipolaron mass about 7 times that of a free carrier 
mass (the bandwidth of a free electron is 8t — 2eV). This 
is not a huge enhancement, since it means that the bipo- 
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FIG. 8: Ground-state bipolaron (a) binding energy, and (b) 
effective mass in units of the free carrier mass vs. Q, for 
various polarizabilities. The full lines correspond to t' = 0, 
dashed lines to t' = -t/2. Here Uh = lOeV. 



laron is about 3.5 times heavier than a single polaron. In 
Fig. 8(b) we plot the bipolaron mass for various a p and 
O values, showing only limited variation over a wide pa- 
rameter range. The higher energy d-pair is heavier, with 
a much narrower bandwidth. 

If we include 2 nd nn hopping t' — —t/2, two effects are 
apparent. First, the two-polaron continuum is pushed to 
higher energies, effectively increasing the binding ener- 
gies of both bound bipolaron states. (The binding en- 
ergies for the GS bipolaron are shown in Fig. 8(a) for 
various parameters). Second, the d-state becomes the 
ground-state. This is not surprising, since the 2 nd nn 
hopping links directly the two nn bipolaron configura- 
tions which give the bulk contribution to these eigen- 
states. Since t' < 0, this mixing raises the energy of an s 
state and lowers that of a d-state. Thus, if the effective 
if between the two nn bipolaron configurations is large 
enough, the d state has to become the ground-state. 

Let us now comment briefly on the triplet eigenstates. 
Since there is no onsite triplet bipolaron configuration, 
their energies are essentially identical to the energies of 
the singlet eigenstates in the limit Uh — ► oo. This im- 
plies that at this level of approximation and for large 
enough Uh, singlet and triplet bound bipolarons are al- 
most degenerate. However, it is well known that second- 
order perturbation theory produces an exchange interac- 
tion which strongly favors the singlet eigenstates (this is 
the interaction responsible for binding the SI Holstein 
bipolaron, see for example Ref. 22). It follows that in 
reality, the ground-state must be a singlet. 

The binding energies shown in Fig. 8(a) are substan- 
tial, even without inclusion of this singlet exchange en- 
ergy. As discussed in the supplementary material, re- 
laxation of our approximations (that only As ions nn 
to a charge are polarized and that dipolc-dipolc inter- 
actions are ignored) further enhances U\, and therefore 
these binding energies, to several eV. Yet more enhance- 
ment is expected if we include even higher orbit als than 



5s when describing the As ion polarization. All this sug- 
gests the appearance of pre- formed pairs well above room 
temperature. On the other hand, in our model Hamilto- 
nian we have ignored a nn repulsion energy which comes 
from the bare Coulomb interaction itself (in reality, this 
nn Coulomb repulsion will be reduced by other screening 
mechanisms, such as the bond polarizabilities involving 
the As Ap and Fe 3s and Ap states). This nn repulsion will 
decrease U\ substantially, and may even change its sign, 
making nn bipolaron pairs unstable. It is difficult to ob- 
tain accurate quantitative estimates of all these terms, to 
find out whether bound bipolaron states exist and what 
are their binding energies. 

It is important to point out that the presence of this 
As-mediated nn attractive interaction U\ may be essen- 
tial even if preformed bipolaron pairs do not exist. This 
would put these materials in a BCS-like framework, with 
the phonon glue replaced by a virtual excitonic glue. Be- 
cause Ui is so large one does not really need much of 
a retardation effect to overscreen the bare Coulomb re- 
pulsion. In favor of this scenario is experimental data 
indicating higher T c in samples with shorter Fe-As dis- 
tances and smaller Fe-As-Fe angles, which is precisely 
what increases the value of U\ (more speculation on this 
is presented in the supplementary material). 

However we do not want to rule out the presence of pre- 
formed pairs of the kind discussed above. There is ample 
evidence that singlet pairs could exist to quite high tem- 
peratures in these systems. For example, the magnetic 
susceptibility [23] shows the same strong increase with 
temperature above the spin density wave (SDW) or su- 
perconducting transition temperature. This is very diffi- 
cult to explain if we start from a free single-particle-like 
picture, but easy to understand if we assume the presence 
of singlet preformed pairs well above T c , with a binding 
energy of the order of lOOmeV or more. The magnetic 
susceptibility of such a system would increase with in- 
creasing T with an activated kind of behavior. NMR 
Knight shift data also displays this kind of behavior [24] . 

We therefore suggest to take the scenario of preformed 
singlet bipolarons seriously. We note that a supercon- 
ducting state would then behave more like a Bose Ein- 
stein condensate. At first glance one might think that its 
T c would have to be very low, however we note that these 
bipolarons are very light, with a mass which is about 3-4 
times the single polaron mass, and thus very high con- 
densation temperatures could be expected. The other 
rather interesting aspect to consider is that perhaps this 
scenario might also explain the low amplitude SDW ob- 
served at low dopings, as being a different kind of conden- 
sate of singlet bipolarons due to rather strong exchange 
interactions. The low amplitude would be a result of 
the pair-wise singlet formation tendency competing with 
the pair-pair exchange interactions which would favor the 
SDW. We propose to study these issues next. 
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I. THE MODEL: LATTICE AND HAMILTON! AN 

The idealized FeAs layer of interest is sketched in Fig. 1, which shows the cubic lattice 
(lattice constant a = 2. 8 A) which has the Fe in the center of the cubes, on a simple square 
lattice. The As are divided equally between the top and bottom layers, as indicated. Each 
Fe has 4 nearest neighbor As (at a distance R = a\/3/2 = 2AA), arranged in a perfect 
tetrahedron, two in the upper and two in the lower layer. (The real structure has somewhat 
distorted tetrahedra). 

In our model, we assume that the valence electrons occupy 3d orbitals of the Fe atoms. We 
are interested in excess charges on the Fe sites, be they electrons or holes, and we describe 
them in an effective single band model: we assume that they can hop between nearest 
neighbor Fe, with a hopping integral t. We will also consider separately the effect of second 
nearest neighbor hopping, with a hopping integral t'. There is very little hybridization of 
these 3d orbitals with As orbitals, and we ignore it altogether [1] . On-site Hubbard repulsion 
for these charges, characterized by an energy Uh, is also included. The value of t,t', Uh and 
of the other parameters introduced below is discussed in the next subsection. Thus, the 
hopping Hamiltonian for the extra charges on the Fe sites is: 

H Fe =f + f' + U=-t J2 ( C l,cr C .J> + h - a ) ( C i,a C j,a + h - C ) + U "*T™U ( l ) 

where c| j(J creates an extra charge on the Fe site i with spin a, h ia = c\ a c ia and 

denote summation over the nearest, respectively second nearest neighbor Fe sites. A more 

detailed model taking in consideration the various relevant 3d orbitals is of course possible, 




FIG. 1: 3D sketch of the FeAs layer. The lattice constant is a and the angle 9 between the Fe-As 
direction and the z-axis is indicated. 
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but not necessary for our purposes. We avoid such complications in order to make more 
transparent the essential new physics proposed in our model. 

Ab-initio calculations find that the anions are effectively in the As 3 " state, with fully 
occupied 4p bands lying well below the Fermi energy [1]. In reality, these As 4p states are 
strongly hybridized with Fe 4s and 4p orbitals (this hybridization is responsible for the fact 
that the effective As charge inside the muffin tin is much less than the ionic charge of —3, 
[1]). Since this mixing is strongly dependent on the Fe 3d occupation, it will provide a 
screening mechanism for the bare Coulomb repulsion, in particular for the value of Uh- In 
the following we consider that such screening is already included in Uh- For simplicity, we 
also refer only to the full As 4p orbitals, ignoring their hybridization with the Fe 4s and 4p 
orbitals (this could also be accounted for in a multi-band model, but is of limited relevance 
and we ignore it for clarity). 

In our model, the only role of these As ions is that they become polarized due to the 
electric field created when an extra charge (electron or hole) is on a nearby Fe atom. This is 
due to transitions of As electrons from the filled 4p shell into (primarily) the empty 5s shell, 
and results into a dipole moment (polarization cloud) pointing towards the Fe atom which 
hosts the extra charge. For simplicity, we only consider the polarization of the 4 As atoms 
closest to the Fe hosting a charge and we ignore the dipole-dipole interactions between these 
As polarization clouds. The quantitative importance of these approximations is discussed 
towards the end of this supplementary material. 

Three more observations are in order. First, it is more convenient to describe the As 
polarization clouds not in terms of an electron being excited from the filled 4p orbitals into 
the empty 5s orbital but, equivalently, in terms of a hole being excited from the 5s into the 
4p orbitals. Thus, each As will be described by 4 operators: s^pl ,p], , which create a 
hole of spin a in the respective orbital. The ground-state, then, is s|sj|0) and the p-orbitals 
lie at an energy Q above it. If we measure energies above the s-orbital, the Hamiltonian of 
the unperturbed As ions is: 

n As = n £ p\ Xa Pi,x,„ (2) 

i,\,cr 

where i is the location of the As atom (see below), X — x,y,z and a is the spin of the hole. 

The second observation is that this Hamiltonian ignores hopping of charges between As 
atoms. This is a reasonable approximation given the large distance between neighbor As 
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atoms, which results in rather narrow As bands, as demonstrated by ab-initio calculations. 
This is why we model the As within this rather atomistic picture. 

Finally, the lattice illustrated in Fig. 1 suggests that we should consider a unit cell with 
(a minimum of) 2 Fe atoms and 2 As atoms, one in the top and one in the bottom layer. 
However, given our assumptions, in our model it makes no difference whether half the As 
are above and half are below the Fe layer, or whether they are all in the same plane (for 
example, in the layer below the plane of Fe). Mathematically, this follows because for the 
As in the layer above, one can redefine the p z operators as p z — > — p z , which is equivalent 
to a local flipping of the direction of the z-axis. This is very convenient, since it allows us 
to use a simple unit cell with one Fe and one As in the basis, as illustrated schematically in 
Fig. 2. However, it is important to point out that if one wants to include, for example, As 
dipole-dipole interactions, than one needs to carefully consider the true orientations of the 
As polarization clouds with respect to one another. 

Figure 2 also indicates our choice for indexing the lattice sites, with Fe; and Asj being the 
unit cell of the simple square lattice. If the charge is on Fe^, the polarization of Asj and the 
other 3 As neighbors is pointing in its direction. It is therefore more convenient to use the 
1,2 axes orientation for the in-plane directions (see Fig. 2), and 3 = z for the out-of-plane 
direction, in terms of which p x = ^{p x — P y ),P2 = ^75 (Px + P y ),P3 = Pz and Has remains 
formally identical to that of Eq. (2), except that now A = 1, 2, 3. 

With these definitions, the cloud of Asj is pointing towards Fe, in the direction — sin 6e 2 + 
cos^, where cos# = ^ comes from simple geometry illustrated in Fig. 1. Of course, Fe« 
also interacts with Asj_ z , Asi- x - y and Asi- y . Their polarization vectors are obtained by 
appropriate rotations from the one just listed. 

The interaction Hamiltonian is, then: 

ftint = Q^hi \s\ a (- sin 9p ij2 ,a + cos 9pi^ a ) + s\_ xa (sin Qp^ xXa + cos 9pi- x , 3 ,a) 

i,cr 

+ 4-x-y,a ( Sin 0Pi-x-y,2,a + COS 6pi- x -y,3,a) + s}-i/,«t (- SU1 Opi-y^o + COS 0Pi-y, 3 ,a) + h.C. 

(3) 

where the sum is over all unit cells i in the lattice. This interaction describes the s — > p 
transitions and resulting polarization of the As clouds (with the appropriate orientations) if 
extra charges, counted by hi = + n^, are on Fe^ site. 
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FIG. 2: Top view of the Fe-As structure. Several sites are indexed. The orientations of the in-plane 
axes 1,2 and x, y is also indicated. The Fe and As layers have different depths in the z direction. 

Our model Hamiltonian Ti = Hf c + Has + T~imt is thus characterized by 5 parameters, 
t, if , Uh, fi, 9- We now discuss the choice of their values. 



A. Parameters 



Ab-initio calculations find the Fe 3c? dispersional part of the bandwidth to be roughly 
W = 2eV [1]. Our 2D hopping model has a bandwidth W = 8t irrespective of the value 
of t', thus t = 0.25eV is a typical value. For t', we will consider both the case if — and 
t' = —t/2 = — 0.125eV, which is the appropriate value for hopping between d xy and the like 
orbitals which contribute most near the Fermi level. The Hubbard repulsion Uh is certainly 
a large energy. We will use it as a parameter, and show that our results are very weakly 
sensitive to its precise value once Uh > 8 eV or so, which is certain to be the case. 

The ab-initio calculations also suggest a difference from the top of the As 4p bands to the 
bottom of the As 5s bands of about 4eV. This value provides a lower limit for the split Q 
between the s and p orbitals. A more appropriate measure would be the distance between 
the mid-points of these bands, but this is harder to estimate. We use Q = 6eV as a typical 
value, and also study the effects of larger Q values on the results. This is also reasonable 
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because higher-energy orbitals of the As will be polarized as well, so one can think of these 
s and p orbitals as modeling effectively the polarization of the whole atom. 

The final parameter needed is g. Since this characterizes the s — p hybridization when 
an As polarization cloud is created, we should extract it from the As polarizability a p . Its 
value is not known precisely, but the polarizability is typically equal to the volume of the 
atom. Since As is a big atom, its polarizability is estimated to be a p = 10 — 12 A 3 . 

Consider the polarization cloud formed on an As ion, due to a charge q placed at a 
distance R away, in the direction characterized by the versor e. The electric field at the 

— * 

As site is then E — j^e. Note that we use an unscreened value for this electric field. The 
reason is that we are interested in behavior related to the high-frequency dynamic electronic 
polarizability, and not that due to slow vibrational motion of the As themselves. As we show 
below, the effective masses of the polarons and bipolarons are fairly small, meaning that they 
are highly mobile. The As ions are fast to acquire a polarization cloud, but since they are 
very heavy, they simply cannot move fast enough to follow the fast polaron dynamic. This 
separation of energy scales allows us to ignore lattice dynamics effects in the following. The 
interaction between the As and the electric field is TCmt = _ V ' E, where p is the As dipole 
moment operator. In the second quantization, and if we restrict ourselves only to s and p 
orbitals to describe the As, its dipole moment is 

P = J2( s \ ef \P^) (4pa, ct + /i.e.) (4) 

A, (7 

where X — x,y,z or A = 1,2,3 (either basis is equivalent), and we already took in consider- 
ation the fact that there are no matrix elements of the dipole moment between orbitals of 
the same level ((s|r|s) = 0, etc). Moreover, (s\er\p\) = ee\dAs, since the operator r has a 
non-vanishing matrix element only in the direction e\ of the orbital p\ involved. Thus, a As 
is a characteristic "size" of the ion, defined as: 

a As = (s\x\p x ) = J drcf)* x {f)x(f) Px {f) = (s\y\p y ) = (s\z\p z ) (5) 

In terms of this, we find: 

ftint = e* s lP\° + h - c ) ( 6 ) 

A,cr 

This is exactly the form used in Eq. (3), except there we have explicitly written the various 
projections e ■ e\ in terms of the angle 9 appropriate for our geometry. It follows that 
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g = —edAsq/R 2 - As we show in the following, the energetics only depends on g 2 , i.e. it 
is the same for both electron or hole extra charges on the Fe atom. Since the sign of g is 
irrelevant, from now on we assume that the extra charges are electrons, q = — e, so that: 

9=^>0. (7) 

This shows that g is an energy equal to the typical induced dipole eaAs times the applied 
electric field e/R 2 . 

For simplicity, in the rest of this discussion let us denote p a = I] A e • e\p\ j(T , i.e. it is 
the annihilation operator for the p orbital oriented in the direction e of the applied electric 
field. Of course, one can form two other linear combinations of the p x ,Py,Pz orbitals which 
are orthogonal to this direction. These two do not couple to the electric field and do not 
participate in the polarization cloud, thus their energetics is trivial. The interesting part for 
the problem of the As cloud in the presence of the extra charge is therefore described by the 
Hamiltonian: 

h = nY / plp. + gY,( s lp* + h - c -) ( 8 ) 

This Hamiltonian is trivial to diagonalize. Its ground state has the energy 

E cloud = n- ^ 2 + 4g 2 (9) 



and the eigenstate 



where 



and 



|GS)= 7 }7l|0> (10) 
7^ = cos as\ — sin ap\ (11) 



cos a 



sm a 



i/ n 



\2\ Vfi 2 + 4g 2 , 



^-TFfv) (12) 



Indeed, note that if g — this is reduced to the As ground state (holes in the s orbitals, 
a = 0°) whereas for g — > oo, the 7 operators describe the maximally polarized linear 
combination (s — p x )/y/2 and a = 45°. 
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The second eigenstate of h has an energy E exc = Q + ^Vl 2 + Ag 2 and corresponds to 
filling up the orbitals r)l = sinasj. + cosctpj., i.e. of creating a polarization cloud oriented 
antiparallel to the electric field. As already discussed, one can also excite the electron in the 
p orbitals orthogonal to the direction of the electric field, the energy of those states being 
simply Q. Given that Q is a large energy, all these excited states are well above the ground 
state and will be ignored in the rest of the calculation. 

Finally, the dipole moment induced on the As atom is: 

(p) = (GS\e-p\GS) = 4ea As sin a cos a = J^2 + \ g 2 ( 13 ) 

Remember that Eq. (7) revealed that g is proportional to the applied electric field E = e/R 2 . 
This shows that in the limit g <C Q we can approximate 

{p} ~ ~Jr~ ~ ~oT E 

and thus identify the polarizability from this linear regime where (p) = a p E. This allows us 
to extract a^ s in terms of a p and then, through Eq. (7), to calculate: 

n Qe 2 

Using a p = 10A 3 and £1 = 6eV gives the estimate g = 2.5eV and thus g/Q = 0.4, suggesting 
that non-linear effects described by our model are starting to become important, however 
they are not dominant yet. For example, with these numbers, the energy of the polarized 
As cloud from Eq. (9) is -E c i ou d = —1-78 eV, whereas in the linear regime §/0 <C 1 we would 
have -E c ioud ~ — ^ = — 2.05eV. As a consistency check, also note that the linear expression 
for the energy of the cloud can be rewritten as -Ecioud ~ — ^ = —^ E ^, as expected to be 
the case in this regime. 

To conclude, we will use t = 0.25eV, t' = or -0.125eV, Q = 6eV and g = 2.5eV as 
typical numbers. As discussed previously, we will also allow Q to vary and we will then 
choose g according to Eq. (14), so that the polarizability stays constant. The corresponding 
g/Q ratio is shown in Fig. 3: for a wide range of values of Q, g/Q is roughly 0.4-0.5. 



II. THE SINGLE ELECTRONIC POLARON 



We would like to find the eigenstates of the total Hamiltonian H = Hf c + 'Has + 'Hint in 
the case when there is a single extra charge at one of the Fe sites. Of course, this charge 
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FIG. 3: Ratio g/Q vs. 17, for a fixed polarizability a 



p 



10 A 3 (full line) and 12l 3 (dashed line). 



The dot indicates our typical choice of parameters. 

can hop around and will also polarize the As atoms in its vicinity. We call this composite 
object (the charge dressed by the surrounding polarization clouds) an electronic polaron. 

To find its low-energy spectrum, we will use perturbation theory in the hopping Hamil- 
tonian. The reason this is a valid approach is that t is by far the smallest energy scale in 
this problem. A more detailed justification is provided below. 



In the absence of hopping the electron will always stay on the same site, say at Fe^, and 
will polarize its 4 neighbor As. Since we ignore dipole-dipole interactions in our model, this 
simply leads to the appearance of 4 polarization clouds like the one discussed in the previous 
section, each pointing towards the Fe atom. This is illustrated in Fig. 3. 

The ground-state energy of the static polaron is then, immediately: 



spin-up or spin-down hole of one of the As atoms polarized in a direction perpendicular to 
its local electric field. These excited states will be ignored in the following, since t <C Q. 

To characterize the polarization clouds, we introduce the following operators [analogs of 
Eq. (11) for the specific geometry here]: 



A. Zero order perturbation: t = 




(15) 



The lowest excited state is at an energy 



g ft + y/Q 2 + 4g 2 above Ep t cs, and has either the 



7j,A,±, CT = cos «4<x-sina 



(± sin 9p[ Xa + cos 9p{ 3ja ) 



(16) 
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FIG. 4: Pictorial illustration of a static electronic polaron. The electron residing on the central 
Fe atom polarizes the holes on the nn As sites into the 4p orbital pointing towards the Fe. More 
distant As are assumed to be in the 5s ground-state. 

where here A = 1,2 shows if the in-plane component of the polarization cloud is oriented 
along axis 1 or axis 2, and the ± sign indicates whether the in-plane polarization is parallel 
or antiparallel to this axis. The angles a and 9 were denned in the previous sections. 
In terms of these, the ground-state of the system with the electron at site Fe^ is: 

1*0 = 410 (i7) 

where \i) describes the state of all As atoms when the electron is at site Fe, (since there 
is a single such electron, its spin is irrelevant and we drop its index here). Following our 
discussion in the previous section, we have 

I*) = nil2-Wi-xX+,v7i-x-y,2,+Wi-yA,-J GS )i ( 18 ) 

(7 

where 

\GS)i= n °u°)> ( i9 ) 

\j-i\>R,a 

i.e. it describes unpolarized (ground-state) As atoms at all sites which are not nearest 
neighbor to the site Fe^ containing the extra charge. Thus, the product of 7 operators 
describes the polarization clouds on the 4 neighbor As sites (both spin-up and spin-down 
holes are polarized), and \GS)i describes the remaining holes on the unpolarized As sites. 
Of course, this ground-state is highly degenerate since the electron can be at any Fe 
site. The hopping t lifts this degeneracy. 
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FIG. 5: Single polaron hopping in the x-direction. The hopping integral is renormalized by the 
overlap between the polarization clouds in the initial and the final position. The unpolarized As 
atoms are in their ground-state s-orbital. 

B. First order perturbation 



In the subspace generated by all the ground-states described in the previous section, 
we can form eigenstates which are invariant to translations in the usual fashion: 



l*s> = E 



N 



(20) 



where the lattice has N x N = N 2 unit cells and we assume periodic boundary conditions. 
As usual, the quasi-momenta k are denned inside the first Brillouin zone of the simple square 
lattice <C k 7 , ^ . 

Since the hopping Hamiltonian is diagonal in this basis, ($g|T tot |$p) oc it follows 
that |$?) are the eigenstates in first-order perturbation theory. Their energies are: 



E P (k) = Ep, GS + ($r|T tot |$n) 



(21) 



with the second term describing the polaron dispersion. Straightforward calculation reveals 
that: 

E P (k) = Ep tGS - 2t eS [cos(k x a) + cos(k y a)} - 4t' e{i cos(k x a) cos(k y a) = E P ^ GS + e eS (k), (22) 

in other words the dispersion is identical to that of a free particle, but with renormalized 
hopping integrals (or mass). The renormalization is due to the overlap of the corresponding 
As clouds as the electron moves from one site to a nearest or second-nearest neighbor site, 
as illustrated schematically in Fig. 5 for the former case. These overlaps result in: 
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t' eS = t'(i\i + x + y) = t' 

The results are plotted (in units of bare t, respectively t') in Fig. 6, vs. Q, for the two 
values of polarizability a p = 10 and 12 A 3 . As expected, in the limit of large f2, which 
implies g/Q — > 0, the ratios increase towards unity because the polarization clouds become 
very small. As Q decreases, the effective coupling g/Q increases and the ratios (and therefore 
the bandwidths) decrease. This is typical polaronic behavior, however note that even if we 
are in the strong coupling limit g/Q — > oo, the effective hopping is lowered to at most 
t eS /t = 1/3 4 = 0.012, t' eS /t' = l/(2 6 3 2 ) = 0.0017. While these are, of course, small values, 
the renormalized hoppings are not becoming exponentially small with increased coupling, 
as is the case for polarons with phonon clouds. The reason is that within this model, the 
polarization clouds saturate to a maximum possible value and therefore the overlap cannot 
become arbitrarily small. 

The effective polaron mass renormalization is inversely proportional to the t e g/t ratio 
(this statement is true if t' = 0. If t' is finite both ratios determine the new effective mass, 
however they are roughly equal). Thus, these polarons remain rather light objects, with a 
mass up to a factor of 3 or less larger than the bare band mass, if Q is not too small. The 
dots in Fig. 6 indicate our chosen values, for which the polaron mass is roughly 2.2 times 
the mass of the free band electron, and £ e ff = 0.45 x t ~ O.lleV. 

While not huge, this renormalization of the effective hoppings/mass is important since it 
further supports our use of the perturbation theory. First order perturbation is reasonable 

0.7 
0.6 
0.5 
0.4 
0.3 

°' 2 4 5 6 7 8 """4 5 6 7 8 
Q (eV) Q (ev) 

FIG. 6: t e g/t (left) and t' cS /t' (right) vs. for different As polarizabilities. The dots indicate the 
main values used here. 
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only assuming that half the effective bandwidth is much smaller than the distance to the 
next set of excited states, which were ignored in this calculation. As already discussed, those 
excited states are at \ f2 + \JVt 2 + Ag 2 « 5eV and higher energies, whereas 4t cff 0.45eV. 
These numbers suggest that first order perturbation theory should be quite accurate for 
parameters in the regime of interest to us (the higher the Q value used, the more accurate the 
perturbation theory is). A second order calculation (which is possible, but rather involved) 
should only bring in relatively small corrections. Moreover, those corrections may well be of 
the same order as other ignored small terms (like dipole-dipole interactions between the As 
clouds) and therefore one would need to reconsider all these approximations carefully before 
going to a higher order in perturbation theory. 

In conclusion, our calculation shows that, to first order, the only effect of the formation 
of the electronic polarons is to scale down the hopping energy scales by a factor of around 
2.5 (with slight differences for nearest vs. 2nd nearest neighbor hopping), but otherwise 
the dispersion is similar to that predicted by LDA for a free charge. ARPES measurements 
indicating precisely such behavior (in LaOFeP) have been reported already in Ref. 2. In 
our theory, ARPES would measure a quasiparticle peak at low energies, with a quasiparticle 
weight 

Z = \(^ k \cl\GS)\ 2 = [cosar (25) 

which is independent of momentum (second and higher order corrections to the wavefunc- 
tions may bring in some small momentum dependence). For our typical values, this implies 
Z = 0.38. The remaining spectral weight should be observed at much high energies typical 
of the excited cloud states, of the order \ Q + ^/Q? + Ag 2 = 5eV or more. 



III. THE BIPOLARON 

Within the same model and with the same approximations (first-order perturbation the- 
ory) we will now calculate the eigenstates of the system when 2 charges are present on the 
Fe sublattice. In particular, we are interested to see whether they become bound into a 
so-called bipolaron, or a state of two free polarons is energetically more favorable. 
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FIG. 7: Polarization clouds for a 2nd nearest neighbor electronic bipolaron. The central "shared 1 
As atom has a polarization different from that of the usual polaron clouds. 



All the eigenstates of the two charges must be either a singlet or a triplet in the spin 
space. Since we expect a singlet to be the ground-state (this expectation is verified by our 
calculations) we focus on calculations for singlet eigenstates. Most of the calculations for 
triplets are very similar, therefore we will only point out where there are major differences 
between the two cases. 

For convenience, we introduce a singlet creation operator, denned as: 



Within zero-order perturbation theory (t = 0) the charges cannot move, being pinned at 
their initial locations. Depending on the distance between them, different types of polariza- 
tion clouds can form with different energies. We will now calculate these energies. 

If the two electrons are 3rd nearest neighbors or further apart, \i — j\ > 2a, then their 
polarization clouds do not overlap (there is no As atom that is a nearest neighbor for both 
of the charged Fe sites). Therefore each charge forms the same 4 clouds as described for 
single polarons, and the energy is twice the energy of a polaron: 



If the electrons are 2nd nearest neighbors, they share one As atom, as shown in Fig. 7, 
so only 7 polarization clouds are formed. Six of these are precisely like the ones discussed 



A. Zero order perturbation: interaction energies 




(26) 




(27) 
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FIG. 8: C/2 vs Q, for two values of the polarizability. Because U2 > everywhere, it follows that 
it is energetically more favorable for electrons to stay at distances \i — j\ > 2 than to be second 
nearest-neighbors. The dot indicates our typical value. 

previously, since each of those As atoms is in the field created by a single charge. The 
cloud on the central As atom, however, is different, since it is created by the sum of the 
electric fields from the two charges. Given the geometry of the system, this new cloud is 
polarized along the z-axis, and this z-axis polarization is twice as big because two charges 
contribute to it. Thus, the only change in calculating the energy and eigenstate of this cloud 
is to replace g — > 2g cos 9 = 2g/y/3 in Eq. (8), and of course the generic p orbital is now 
the p z = P3 orbital. As result, the energy of this cloud is changed to Q — ^JQ 2 + ^9 2 an d 
therefore the total energy of the 2nd nearest neighbor static bipolaron is: 



Ebp,2 — 6 



Q - Jtt 2 + Ag 



(28) 



This energy is different from the energy of two free polarons by an amount 



U 2 = E BP>2 - E BP}OC = - n - + y (29) 

This is plotted in Fig. 8 as a function of Q, for the two values of As polarizability. In both 
cases, it is positive and a slowly varying function of Q. The 2nd nearest neighbor effective 
interaction is repulsive because the new cloud, although bigger than a single original cloud, 
is not bigger than two such clouds. As a result, it is energetically more favorable to have two 
free electronic polarons situated far apart, than to have a second nearest-neighbor bipolaron. 
This can be shown to be true for any value of the polarizability and of Q (i.e any g/Q ratio), 
for this lattice structure. 
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FIG. 9: Polarization clouds for a nearest neighbor electronic bipolaron. The two central "shared" 
As atoms have a polarization different from that of the usual polaron clouds. 

Interestingly, this interaction can be made attractive by changing the lattice structure. 
Focusing on the weak-coupling regime of small ^ and using a linear approximation, the 
energy of the new cloud is Q — \/Q 2 + lQg 2 cos 2 9 pa —8 cos 2 9g 2 /Q, while the energy of two 
of the original clouds is pa —4g 2 /Q. The 2nd nearest neighbor pair becomes energetically 
favorable when 2 cos 2 9 > 1 — > 9 < 45°. This can be achieved by making the lattice 
tetragonal, i.e. by increasing the distance between the Fe and As layers by a factor of \/2 or 
more. Of course, this would also have the effect of decreasing all the polarization energies, 
since the electric fields would be smaller due to the increased Fe-As distances 

The ground-state eigenstate of this new cloud can be obtained from the general formula, 
after performing the substitutions discussed above. This leads to the appearance of a new 
mixing angle: 



cos f3 



sin (3 



\ 



1 



n 2 + f g 2 



\ 



n 



(30) 



in 2 + fg 2 / 

which appears in the definition of the new cloud. For instance, assuming that this new cloud 
appears at Asj site, then its creation operator for a given spin is: 

cos/3st CT -sin/3pt 3jCT (31) 



1 1,(7 



A similar analysis applies for a nearest neighbor bipolaron. In this illustrated in 

Fig. 9, a total of 6 As atoms are polarized. Four of them have clouds of the original type, 
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FIG. 10: Left panel: Effective nearest neighbor bipolaron interaction U\ as a function of f2 for two 
polarizability values. Right panel: U\ vs. g/£l for Q = 4,6 and 8 eV. 

while 2, which are neighbors to both Fe atoms, have yet another type of polarization cloud 
which lies in the plane perpendicular to be bipolaron axis. Straightforward geometry shows 
that this polarization cloud has magnitude larger than an original cloud, therefore one 
has to replace g —>■ gJ\ (or 4g 2 — > y<? 2 ) to find the energy of each new cloud. 
Thus, the total energy of a static nearest neighbor bipolaron is: 



E 



BP,1 



- J 9? + Ag 2 



+ 2 



n- Jn 2 + 32 g 



(32) 



and the difference between this and the energy of two unbound polarons is: 



Z7i = E BP ,i - E BPM = AJQ? + Ag 2 - 2tt - 2\Kl 2 



32 



-r 



(33) 



This effective energy is plotted in Fig. 10 both as a function of Q for a fixed polarizability 
(left panel) and as a function of g/Q (right panel) for several values of Q. 

The first interesting observation is that this effective interaction is attractive, favoring 
a bound nearest neighbor bipolaron over two independent polarons. As shown in the right 
panel, this is only true for small g/fl < 1 values (if g/Q becomes large enough, the clouds 
saturate to their maximum allowed values and the two new clouds are energetically less 
favorable than four original clouds). Remarkably, the strongest binding energy is for g/Q « 
0.45 for reasonable values of Q, which is close to the values we expect for this parameter 
(see Fig. 3). From this point of view, we can say that these lattices are already very close 
to being fully optimized. As discussed for the case of the 2nd nn bipolaron, it is possible 
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to further change this interaction by changing the lattice geometry: a smaller 9 angle will 
again increase this binding energy, but this is balanced by a decrease if the Fe-As distance 
is increased. 

The creation operators for these new clouds are: 

^p\,\,a + Pi,3,a 



li,±,o =cos74, CT -sm7- 



V2 



(34) 



where A = x or y is the direction perpendicular to the bipolaron axis, and the ± sign shows 
whether the cloud is parallel/antiparallel to the A axis. The new mixing angles are: 



cos 7 



\ 




VI 2 + fg 2 



sin 7 = 



1 



1 



^fg 2 



(35) 



The analysis to this point would be identical if the static bipolarons were in a triplet state. 
The last case, namely that of an onsite bipolaron, is however only possible for a singlet state. 
In this case, there are again only 4 polarization clouds, as for a polaron. However, because 
the charge creating the electric field is doubled, we now replace g — > 2g. As a result, the 
energy of the on-site bipolaron is: 



E 



BP,0 



The effective on-site interaction is then: 



Q - ^/Q 2 + 16g 2 



(36) 



Uq = Uh + E B pfl — E B p,oo = Uu — 4 



tt 2 + lQg 2 + - 2JQ 2 + Ag 2 



(37) 



where Uh is the on-site Hubbard repulsion between the two charges. Clearly, the polarization 
of neighboring As atoms always acts to screen out the on-site repulsion, and this effect can 
be very significant. Indeed, in Fig. 11 we plot this renormalization energy, which is several 
eV for typical parameters of interest to us. 

One more comment is in order. In the linear regime, this renormalization energy would 
be — 16g 2 /f2 ps —lQeV, which is similar to the values used in Ref. [1]. The smaller values 
we find here are a direct consequence of the non-linear effects, which are accounted for our 
this model. 



19 



O 

> 




FIG. 11: Renormalization of the on-site interaction Uq — Uh vs. Q, for two values of the polarizabil- 
ity. This is always a large attractive contribution, which significantly screens the on-site Hubbard 
repulsion. The dot shows the typical value we use. 

The cloud operators for the on-site bipolaron are similar to the original polaron cloud 
operators: 



=t 

7i,A,±,<T 



cos 5s 



t 



sin 5 



± sin I 



cos I 



4,3,o- 



(38) 



where A = 1 or 2 depending on the in-plane alignment of the cloud, and ± reflects its 
orientation. The only difference is the new mixing angle: 



cos 5 



sin 5 



\2 



1 



\2 



1 



n 



(39) 



To summarize, this calculation reveals that for the parameters of interest to us, we can 
have a stable on-site bipolaron if the Hubbard repulsion is not too large. The nearest- 
neighbor bipolaron is bound, whereas the 2nd nearest neighbor bipolaron is unstable to 
dissociation into two free polarons. 



B. First order perturbation: effective hopping integrals 

We continue the discussion for singlet bipolarons, and analyze how the hopping Hamil- 
tonian mixes together the low energy bipolaron states sj\ + 5). In this notation, 5 
measures the distance between the two charges, and the state + S) describes the As 
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atoms when there are charges at Fe sites i and i + S. Similar to the states \i) defined for 

— * 

single polarons, the states + 5) have all As atoms which are not nearest neighbors to 
either Fe^ or Fe i+ g in their ground state, while the As atoms in the vicinity of the charged Fe 
atoms are described by the appropriate 7 T ,7 T ,7 or 7 operator, as detailed in the previous 
section. We do not list these expressions explicitly here, as they are rather long and tedious 
but otherwise straightforward. 

From these states, we can define a basis of bipolaron states invariant to translations on 
the N x N unit cells lattice: 

p ik-{Ri + l) x 

M> = E^^i + ,M + 5>. (40) 

i 

Here k is the center-of-mass momentum of the bipolaron pair. It can take the usual equally- 
spaced values inside the Brillouin zone allowed by the periodic boundary conditions. How- 

— * 

ever, some care needs to be taken in the definition of 5, since if we take all its possible values, 
we double-count the states (for example, \k, e x ) = \k, —e x ), etc). This is due to the fact that 
5 is like a nematic vector, with a size and orientation but without a pointing arrow. Let 
5 = 5 x e x + 5 y e y , and for simplicity assume that the lattice size is an odd number, N = 2n+l. 
Then, double-counting is avoided provided that when 5 X = 0, we allow 5 y — 0, 1, . . . , n (be- 
cause of the periodic boundary conditions, the distance between charges cannot be more 
than half the dimension of the system). If 5 X — 1, 2, . . . , n, then 5 y = —n, n. Both posi- 

— * 

tive and negative values are needed in this case, because, for example, the \k,e x + e y ) state 
is distinct from the \k,e x — e y ) state, however these are the only two distinct states of the 
2nd nearest neighbor bipolarons. In conclusion, for a N x iV lattice with N — 2n + 1, there 
are (n+ l)(2n + 1) — n distinct singlet eigenstates corresponding to a given total momentum 

— * 

k (the remaining are triplet states. Similar accounting can be done for lattices with even 
number of sites). 

Because our total Hamiltonian is also invariant to translations, it will not mix states 
with different momenta. Therefore, we can solve the problem separately in each k subspace, 
and only need to calculate the various (k,8'\Tl\k,8) matrix elements. There are three sets 
of such contributions. First, there are "diagonal" contributions, already discussed in the 
previous section. Specifically, (k,5\H\k,6) = (k,5\U + Has + T^int I k, 5) = Ug where U , U\ 

— * 

and U 2 of Eqs. (29), (33) and (37) correspond to 5 indicating, respectively, an on-site, a 

— * 

nearest neighbor or a 2nd nearest neighbor bipolaron state. If \5\ > 2, then Ug = 0, since we 
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measure the interactions energies with respect to the energy of two unbound static polarons. 
The nn hopping Hamiltonian links states with the same momentum but bipolaron dis- 

— * — * — * — * 

tances 5 and 8' which are nearest neighbor, i.e. \8 — S'\ = 1. We therefore also need to 
calculate these matrix elements (k, 8 ± e x / y \H\k, 8) = (k,5± e x / y \T\k, 8). Similarly, for finite 
t' there are matrix elements of T' between 2nd nearest neighbor states with \8 — 8'\ = \[2. 
Taken together, these three sets exhaust all matrix elements which are finite and therefore 
we can diagonalize the resulting matrix and find the energies and eigenstates corresponding 

— * 

to any desired k, within first order perturbation theory. 

— * ^ 

If 8 is large enough, the hopping Hamiltonian T will hop one of the electrons away from 
its clouds, while the second one remains unperturbed. As a result, one expects the matrix 
element to be proportional to the t e s of the single polaron. Indeed, explicit calculation shows 
that if the separation 8 is at least that of a 3rd nearest neighbor bipolaron, then: 

(k, 8 + e x \f\k, 8) = -2t eS cos ^ 

(k, 8 ± e y \f\k, 8) = -2t cS cos ^ (41) 

— * — * 

The cosine factors appear because of the phases in the definitions of the eigenstates \k,S), 

— * 

which change when the singlet separation 8 changes. Given the restrictions on the allowed 

— * 

values of 8, one needs to be a bit careful near the boundary of allowed values. For example, 
based on the definition of Eq. (41), it is straightforward to check that 

\k, S x e x + (n+ l)e y ) = e ikya ^\k, S x e x - ne y ). 

Given the periodic boundary conditions, the allowed values for momentum, if N — 2n + 1, 
are k y a = 2itn y /{2n + 1), and the phase factor linking the two states can only be ±1. This 
sign is important because hopping in the y-direction links the state with S y = n to that with 
8 y — n + 1, however the later one is not amongst the allowed 8 values. The identity above 
links it to an allowed value, but also introduces a phase shift (the sign) which is transferred 
to the corresponding matrix element. A similar situation arises for hopping in the x direction 
when 5 X = n, 5 X = and hopping in the y direction for 5 y = —n. Finally, it is important to 
note that these phase shifts are different for triplet states (they have opposite sign in most 
cases, but not all). Other than this (and the non-existence of an triplet onsite bipolaron, as 
already discussed), essentially everything else is the same for triplet states. 
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Different effective hopping integrals appear when the electrons are closer together. For 
nearest-neighbor hopping, there are four such special cases. The first corresponds to hopping 
from 3rd nearest neighbor to nearest neighbor bipolaron. A straightforward calculation of 
the overlap of the corresponding polarization clouds leads to an effective hopping: 

4 



£3 = t{i, i + e x \i, i + 2e x ) — t cosa(cos a + sin a sin 9) 

/ cos ti sni H \ 1 

X 



cos a cos 7 + sin a sin 7 
The matrix elements for such processes are: 



cos 9 sin 9 
+ — — 



V2 



(42) 



Jv CI 

(k, 2e x \T\k, e x ) = -2£ 3 cos -y- 

— * — * Jv Ob 

(k, 2e y \T\k, e y ) = -2t 3 cos 



(43) 



in other words like the general matrix elements of Eq. (41), except with the appropriate 
value of the effective hopping integral. 

The second special case involves hopping between 2nd and 4th nn bipolarons. In this 
case, the overlap between the corresponding polarization clouds leads to: 

r 1 4 

£4 = t{i, i + e x + e y \i, i + e x + 2e y ) = t cos 6 a cos 2 a + sin 2 a cos 2 9 

x [cos a cos (3 + sin a sin (3 cos 9] 2 (44) 

and the matrix elements: 



(k, 2e x ± e y \T\k, e x ± e y ) = — 2t 4 cos 



(k, e x ± 2e y \T\k, e x ± e y ) = — 2t 4 cos 



(45) 



The third special case involves hopping between 2nd and 1st nearest neighbor bipolarons. 
In this case, the overlap between the polarization clouds corresponding to these two cases 
leads to a renormalized hopping: 



t 2 — t(i, i + e x \i, % + e x + e y ) = t cos 6 a 

( cos 9 sin 9 

x cos a cos 7 + sin a sin 7 
and the matrix elements: 



cos a + sin a cos 9 
2 



n 2 



COS j3 COS 7 + 



sin (3 sin 7 



x/2 



V2 



(46) 



(k, e x ± e y \f\k, e y ) = -2t 2 cos ^ 
(k, e x ± e y \f\k, e x ) = -2t 2 cos 



(47) 
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FIG. 12: Effective hoppings for closely-spaced bipolarons. For these parameters, the ti and t% 
curves are essentially superimposed, therefore we only show tijt. 

Finally, for the hopping between on-site and nearest neighbor bipolaron, the effective 
hopping is found to be: 



t\ = t(i, + e x ) = t cos 4 a cos 4 (a — 5) 
and the matrix elements are: 



cos 8 cos 7 + sin 8 sin 7 



cos 9 sin 9 s 



(48) 



(k, e x \T\k, 0) = — 2y/2ti cos -| 



(k,e y \f \k,0) = -2V2U cos 



(49) 



The extra factor of is because of the on-site singlet normalization, see Eq. (26). Of 
course, for triplet states t\ = 0, since there is no on-site triplet. 

In Fig. 12(a) we plot the values of txjt^jt^jt and t^/t vs. Q for otp = 10A 3 . These 
should be compared with t e s/t shown in Fig. 6. The overall changes are relatively small. 

The situation for second nearest neighbor hopping is similar. If S is large enough, then 
the cloud overlap is the same as for individual polarons and the renormalized hopping is t ; eS , 
giving the matrix elements: 

\k x i ky)a 



(k, S + e x ± e y \T'\k, 5) = —2t' cS cos 



(50) 



Again, special care must be used for values of S near the edges of its area of allowed values, so 
that the periodic boundary conditions are properly accounted for. This is done in a manner 
totally analogous to that explained for nearest neighbor hopping. 
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There are again four special cases of different renormalization of t'. One involves the 
hopping between the two possible nearest-neighbor bipolaron states, which has the matrix 
element: 



(k,e x \T \k,e y ) — 2t 2 



(k x + k v )a, (k x k y )ci 
cos — — h cos — — — 



(51) 



where 

t' 2 = t'(i,i+x\i,i+y) = t' [cos a] 



/ cos 9 sin^" 

cos a cos 7 + sma sin 7 — -\ 

V V2 2 / 



1 4 



cos 2 7 + - sin 2 7 



(52) 

The second involves the hopping between 2 nd and 3 rd nn bipolarons, as well as hopping away 
from a 2 nd nn bipolaron, with the matrix elements: 

{k, 2e x ± 2e y \f'\k, e x ± e y ) = -2f 3 cos i^±^> 



(k, 2e x \T'\k, e x ± e y ) = —2t' 3 cos 
(k, 2e y \T'\k, e x ± e y ) = —2t' 3 cos 



[kx ~F ky)a 
2 

(fca, =F k y )a 



(53) 
(54) 



where t' 3 = t'(i, i + x + y\i, i + 2x + 2y) = t'(i + x + y\i + 2x) is found to be: 

t' 3 = t' [cos a] 10 [cos a cos (3 + sin a sin /? cos 9} 2 cos 2 a + + sin 2 a cos (20) 



(55) 



The third case involves hopping between I s * and 4 th nn bipolarons, with matrix elements: 

(k, 2e x ± e y \f'\k, e x ) = -2f 4 cos i^L±^> 
(fc, e x ± 2e y |T / |A;, e^) = — 2t 4 cos 



{kx i ky)d 



(56) 
(57) 



where £4 = + x\i + 2x + y) equals: 



t' 4 = t' [cos a] 8 cos 2 a + sin 2 a cos(20) 



(cos 9 sin$" 
cos a cos 7 + sin a sin 7 | — 7=- + 



(58) 



' 2 

Finally, the hopping between an on-site and a 2nd nearest neighbor bipolaron (again, 
possible only for singlet states) leads to: 



(jfe, e x ± e y \f'\k, 0) = -2V2t[ cos {h ± ky)a 

2 



(59) 
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where 

t\ = t'(i,i\i, i + x + y) —t' [cos a cos(o; — <5)] 6 [cos(/3 — 5)] 2 . (60) 

The renormalized values t^t^t^ and t' 4 are shown, in units of t', in Fig. 12(b). 

The ensemble of all these matrix elements define the matrix that needs to be diagonalized 

— * 

in order to find eigenstates and eigenvalues for a given total momentum k. The results are 
shown in the main paper. 

IV. ESTIMATES FOR THE ACCURACY OF VARIOUS APPROXIMATIONS 

In the previous sections we argued that first order perturbation theory is reasonably 
accurate for our Hamiltonian, given the values of the various parameters relevant for the 
Fe-based superconductors. However, the Hamiltonian itself already embodies two approxi- 
mations, namely (i) that only As atoms nn to a doping charge become polarized, and (ii) 
that dipole-dipole interactions between the various polarized As atoms are ignored. 

We provide here rough estimates for the accuracy of these approximations. We begin by 
discussing the single polaron. 

Assume that interactions with the 8 As atoms which are 2 nd nn to a charge are also 
included. This supplementary interaction would be described by a Hamiltonian similar to 
Hint of Eq. (3), except it would have a new energy scale g' and new angles describing the 
orientations of these Fe-As bonds. 

From the definition of the interaction energy, Eq. (7), we have: 

g R' 2 

where R' = vTla/2 is the distance from an Fe to a 2 nd nn As atom. This comes about 
because the only difference is in the electric fields created at the two sites, and these are 
inversely proportional to the square of the distances. 

It then follows immediately that the contribution to the static polaron energy of these 8 
clouds would be: 

E corr = 8 (Q - ^ 2 + 4<?' 2 ) , (61) 

which for our typical parameters is E COTT = — 1.2 eV. Of course, considering the polarization 
of even more distant As atoms would lower this energy even more. 
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However, this decrease is more than compensated for by consideration of dipole-dipole 
interactions. To estimate these energies, let us consider only the 6 pair interactions between 
the 4 largest dipole moments on the nn As sites. Consider a pair of two such dipoles, 
with moments pi and p 2 . Their magnitude is given by Eq. (13) and their orientations are 
straightforward to determine, as is the distance d between them. Using these values, we find 
for any such pair that: 



where the numerical factor in front is due to geometrical considerations. For our typical 
parameters, £d-d = 0.66eV and therefore the total for the 6 pair interactions is ~ 4 eV. 
Inclusion of contributions from interactions with the dipoles on the 2 nd nn As atoms will 
decrease this, because nn and 2 nd nn dipoles closest to each other are now roughly parallel, 
not anti-parallel, to each other. 

Taken together and without further corrections, these two energies would renormalize 
Eb,ps to about half the value predicted in their absence. If we keep adding further rows 
of neighbors into the calculation, the polaron energy saturates (slowly) to around -1.6eV 
(of course, if the fields decrease as l/R 2 , the energy slowly become more negative and will 
eventually go to — oo. However one has to take into account that there are other screening 
mechanism that further screen the interaction). Having extended clouds would also lower 
the effective hoppings somewhat more, although we expect this to be a much smaller effect: 
clouds which are far from the electronic charge are changed much if the charge moves by 
one site. As a result, their contribution to the total overlap should be very close to 1. 

Although the renormalization of the static part of the polaron energy is quite substantial, 
these corrections do not have nearly as large an effect on the bipolaron interaction energies 
U , Ui and U 2 - The reason is that these are differences between the true bipolaron energy and 
the energy of two free polarons, and therefore many contributions from both the extended 
clouds and the dipole-dipole interactions cancel each other out. 

For example, consider the effects of these corrections on U\, which is the most interesting 
energy. For two polarons at an infinite distance from each other, the dipole-dipole correction 
is 12 .Ed-d and there are a total of 16 2 nd nn extra clouds. For a nn bipolaron, the dipole- 
dipole energy is roughly 11 -Ed-d, because the two As tetrahedra share one side (of course, 
the two central As atoms have somewhat different dipole moments pointing in somewhat 




(62) 
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different directions, and also there are smaller contributions from interactions between the 
non-central As atoms belonging to different Fe, but 11 E d _ d should be a reasonable guess). 
This configuration would only have 12 2 nd nn extra clouds. Subtracting these corrections 
leads to an overall correction of [q roughly equal to —\E COTX — E d _ d 0. In reality, the 
energy is lowered because the close location of the two charges implies larger electric fields 
at all sites for the nn bipolaron, therefore larger clouds. We find that as we increase the 
number of As neighbors included in the cloud, Lq is lowered from « -1 eV (the value we 
find when the approximations are made) to —2.87 eV if we include 2nd nn and dipole-dipole 
interactions, to —4eV if more and more neighbors are added. Again, one expects that the 
electric fields will be screened over some finite distance which will decide the precise value. 
The point, though, is that relaxing the approximations increases this attraction energy, in 
other words the physics we discuss here becomes more relevant, not less. 

V. LINK BETWEEN BCS-LIKE THEORY AND EXPERIMENTAL DATA 

There have been many works in the recent experimental literature pointing to correlations 
between the critical temperature T c and shorter Fe-As distances and smaller Fe-As-Fe angles, 
see for example Ref. 3. Such correlations would come as a natural consequence in our theory, 
where the "glue" leading to pair formation is the polarization of nearby As ions, and thus 
it is very sensitive to the precise location of the As with respect to Fe. 

As discussed towards the end of our main paper, we currently envision two possible 
scenarios explaining the high-T c in the materials, namely a BCS-like and a BEC-like scenario. 
The BEC-like scenario assumes that the nearest-neighbor attraction Lq is larger than the 
nearest-neighbor bare Coulomb repulsion (not included in our simulations) and that pre- 
formed bosonic pairs appear well above T c . As discussed, there is experimental evidence 
consistent with this scenario. The BCS-like scenario assumes that the nn Coulomb repulsion 
is large enough to overcome Lq, and therefore there are no pre-formed pairs in the system. 
However, like in regular superconductors, the existence of an attraction mechanism (given 
by Lq), when combined with the existence of a Fermi sea and proper inclusion of retardation 
effects, should suffice to lead to Cooper-pair formation at T c . 

Currently, our calculations are not detailed enough to allow us to distinguish which 
scenario is the more likely one. In this section we first speculate on the possible fit of the 
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experimental data to the BCS-like scenario, and then comment on the relevance of the same 
data for the BEC scenario. 

The well-known BCS formula links T c to the typical scale of the attractive potential V in 
an exponential dependence, T c ~ exp (— mn,\ where p is the density of states at the Fermi 
energy. On a logarithmic scale, this implies lnT c = const. — if we assume that the 
prefactor is roughly constant. For simplicity, let us also assume that V oc U\ oc a p E 1 ■ E 2 oc 
a p ^r, where R is the Fe-As distance and 7 is the nn Fe-As-Fe angle. As explained in the 
main paper, this expression gives the linear approximation for the nn (/1 attraction [1], and 
we use it here for simplicity - since we will ignore changes in the density of states for the 
various materials, taking into consideration the small non-linear effects is superfluous. 

In Fig. 13 we show a plot of T c vs. — for a multitude of materials, which have 
polarizabilities, R and 7 values varying over fairly large ranges [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
13, 14, 15, 16]. While there is some spread to the data, it is not inconsistent with a linear 
dependence like that predicted by a BCS formula. Of course, a lot more theoretical work 
needs to be done before one could justify this plot as more than just an accident. 

This sensitivity of T c to a p , R and 7 is not inconsistent with a BEC-like scenario, either. 
Assuming that U\ is large enough that pre-formed pairs exist up to very high temperatures, 
then T c would be related to the on-set of coherence between these bosons. Generally speak- 
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ing, this is achieved when the de Broglie wavelength becomes comparable to the average 
spacing between bosons. Whatever the specific relationship is, it will involve the effective 
mass of these pre-formed pairs. This effective mass depends in a non-trivial fashion on a p , R 
and 7, as discussed in the main paper. This dependence is not only through Ui, but to a 
large degree also through the cloud overlaps determining the various t e g values. Since all 
these quantities are very sensitive to a p , R and 7, this scenario would also be compatible 
with a T c which depends on these quantities. Trying to fit it would, at this point, be even 
more speculative than the fit shown in Fig. 13. 

In any event, the underlying idea is that if the As ions played a small, secondary-type 
role in determining the behavior of these materials, this dependence would be very puzzling. 
In our model, the As anions are the main enablers of the pairing mechanics, therefore this 
sensitivity of T c to specific details of the lattice structure follows quite naturally. 
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